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Distributed Stochastic Market Clearing with 
High-Penetration Wind Power 

Yu Zhang, Student Member, IEEE, and Georgios B. Giannakis, Eellow, IEEE 


Abstract —Integrating renewable energy into the modern power 
grid requires risk-cognizant dispatch of resources to account 
for the stochastic availability of renewables. Toward this goal, 
day-ahead stochastic market clearing with high-penetration wind 
energy is pursued in this paper based on the DC optimal power 
flow (OFF). The objective is to minimize the social cost which 
consists of conventional generation costs, end-user disutility, 
as well as a risk measure of the system re-dispatching cost. 
Capitalizing on the conditional value-at-risk (CVaR), the novel 
model is able to mitigate the potentially high risk of the recourse 
actions to compensate wind forecast errors. The resulting convex 
optimization task is tackled via a distribution-free sample average 
based approximation to bypass the prohibitively complex high¬ 
dimensional Integration. Furthermore, to cope with possibly 
large-scale dispatchable loads, a fast distributed solver is devel¬ 
oped with guaranteed convergence using the alternating direction 
method of multipliers (ADMM). Numerical results tested on 
a modified benchmark system are reported to corroborate the 
merits of the novel framework and proposed approaches. 

Index Terms —ADMM, conditional value-at-risk, demand re¬ 
sponse aggregator, market clearing, stochastic optimization, wind 
power. 


Nomenclature 

A. Indices, numbers, and sets 

T, T Number and set of scheduling periods. 

Nh, Ni Number of buses and lines. 

Ng, Mg Number and set of conventional generators. 

Ma, Ma Number and set of aggregators. 

Mw, Mw Number and set of wind farms. 

Ns, Ms Number and set of wind power generation 

samples. 

Ej Set of end users served by aggregator j. 

Srj Set of smart appliances of residential user r 

served by aggregator j. 

Vjrs Set of operational constraints of appliance s 

of residential user r served by aggregator j. 
T^^s scheduling periods of appliance s. 

k ADMM iteration index. 


Manuscript received March 12, 2014; revised August 21, 2014, December 
20, 2014, and Mai'ch 29, 2015; accepted April 01, 2015. This work was 
supported by the Initiative for Renewable Energy & the Environment (TREE) 
grant RL-0010-13, University of Minnesota, and NSF grants CCF-1423316 
and CCF-1442686. Part of this work was presented at the 2014 IEEE 53rd 
Annual Conference on Decision and Control, Los Angeles, CA, December 
15-17, 2014. 

The authors are with the Department of Electrical and Computer 
Engineering and the Digital Technology Center, University of Min¬ 
nesota, 200 Union Street SE, Minneapolis, MN 55455, USA. E-mails: 
{zhanl22 0 , georgios}@umn.edu. 
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Minimum and maximum power output of 
conventional generator i. 

Ramp-up and ramp-down limits of conven¬ 
tional generator i. 

Fixed base load power demand in slot t. 
Maximum power provided by demand re¬ 
sponse aggregator j. 

Minimum and maximum power consump¬ 
tion of appliance s. 

Total energy consumption of appliance s. 
Start and end times of appliance s. 
Minimum and maximum power flow limits. 
Maximum committed wind power. 
Branch-node incidence matrix. 

Incidence matrices of conventional genera¬ 
tors, wind farms, and DR aggregators. 
Nodal susceptance matrix. 

Matrix relating bus angles to branch power 
flows. 

Branch susceptance matrix. 

Susceptance of line i. 

Vector collecting selling prices in slot t. 
Vector of purchase prices in slot t. 
Tolerance of the ADMM termination crite¬ 
rion using primal feasibility. 

Weight of augmented Lagrangian. 

Weight of CVaR-based transaction cost. 
CVaR probability level. 


C. Decision variables 


pt 

Output of conventional generator i in slot t. 

^jrs 

Consumption of appliance s in slot t. 

pt 

^DRAj 

Total power consumption of aggregator j in 
slot t. 

pt 

Power committed by wind farm m in slot t. 

■n 

A variable in the CVaR-based transaction 

cost. 

6/* 

Vector of nodal voltage phases in slot t. 


Vector collecting Pq for all i € Mg. 

Pdra 

Vector collecting J ^ -Yi. 

Pw 

Vector collecting for all m € Mw 

Pjrs 

Vector collecting for all t G T. 

Po 

Vector collecting 77 and p^, pIu’ 

for all t gT. 

Pj 

Vector collecting pjrs for all r and s. 
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D. Uncertain quantities 

Actual power output of wind farm m in slot t. 
w* Vector collecting for all m G Afw 

E. Functions 

a(-) 

Ujrs{‘) 

m-) 

Fp{-) 

Lp{-) 

nG„, nDRA„, 


F. Abbreviations 


ADMM 

Alternating direction method of multipliers 

CVaR 

Conditional value-at-risk. 

DA 

Day-ahead. 

DSM 

Demand side management. 

DR 

Demand response. 

ED 

Economic dispatch. 

ISO 

Independent system operator. 

LMPs 

Locational marginal prices. 

LOLP 

Loss-of-load probability. 

MC 

Market clearing. 

OPE 

Optimal power flow. 

RES 

Renewable energy sources. 

RT 

Real-time. 

SCED 

Security-constrained economic dispatch. 

SCUC 

Security-constrained unit commitment. 

SAA 

Sample average approximation. 

UC 

Unit commitment. 

VaR 

Value-at-risk. 

WPPs 

Wind power producers. 


1. Introduction 

The future smart grid is an automated electric power grid 
that capitalizes on modern optimization, monitoring, com¬ 
munication, and control technologies to improve efficiency, 
sustainability, and reliability of generation, transmission, dis¬ 
tribution, and consumption of electric energy. Limited supply 
and environmental impact of conventional power generation 
compel industry to aggressively utilize the clean renewable 
energy sources (RES), such as wind, sunlight, biomass, and 
geothermal heat, because of their eco-friendly and price- 
competitive advantages. Growing at an annual rate of 20%, 
wind power generation already boasted a worldwide installed 
capacity of 318 GW by the end of 2013, and is widely 
embraced throughout the world im. Recently, both the U.S. 
Department of Energy (DoE) and the European Union (EU) 
proposed ambitious blueprints towards a low-carbon economy 
by meeting 20% of the electricity consumption with renew¬ 
ables by 2030 and 2020, respectively ||2l, E). 

Towards the goal of boosting the penetration of RES, robust 
and stochastic planning, operation, and energy management 


with renewables have been extensively investigated recently. 
A key challenge of the associated power dispatch tasks is 
to account for the intrinsically random and non-dispatchable 
nature of RES so that total power demand can be satisfied 
by total power supply, while the social cost is minimized. 
Being resilient to communication outages and malicious cyber¬ 
attacks, efficient decentralized algorithms deployed over the 
interdependent power entities are indispensable as well. 

Limiting the loss-of-load probability (LOLP), risk-aware 
energy management approaches including economic dispatch 
(ED), unit commitment (UC), and optimal power flow (OPE) 
were formulated as chance-constrained optimization problems 
in a-iii. Leveraging scenario sampling, a general non- 
convex chance-constrained program can be relaxed and solved 
efficiently as a convex one, which however turns out to be too 
conservative in certain scenarios 0. As an alternative, risk- 
limiting dispatch has been formulated as a multi-stage stochas¬ 
tic control problem il; see also m, where direct coupling 
of the uncertain energy supply with deferrable demand was 
accounted for using stochastic dynamic programming. 

Additional early works relied on the so-termed commit¬ 
ted renewable energy. ED penalizing (under-) over-estimation 
of wind power was investigated in Ha. Worst-case robust 
distributed ED with demand side management (DSM) was 
proposed for grid-connected microgrids d. However, the 
worst-case scenario is unlikely to come up in real-time (RT) 
operations. Multi-period ED with spatio-temporal wind fore¬ 
casts was pursued in ifTSll . The obtained optimal operating 
point though can be very sensitive to the forecast accuracy. 

Turning attention to power system economics, market clear¬ 
ing (MC) is one of the most important routines for a power 
market, which relies on security-constrained UC or OPE 
Independent system operators (ISO) collect generation bids 
and consumption offers from the day-ahead (DA) electricity 
market. The MC process is then implemented to determine the 
market-clearing prices lfT4ll . Deterministic MC without RES 
has been extensively studied; see e.g., Ha-El. Optimal wind 
power trading or contract offerings have been investigated 
from the perspective of wind power producers (WPPs) m- 
MC under uncertain power generation was recently 
pursued as well. As uncertainty of wind power is revealed 
on a continuous basis, ISOs are prompted to undertake cor¬ 
rective measures from the very beginning of the scheduling 
horizon ll22l . One approach for an ISO to control the emerging 
risk is through the deployment of reserves following the 
contingencies 12^ . Electricity pricing and power generation 
scheduling with uncertainties were accomplished via stochas¬ 
tic programming 12^ . ||25l. In addition, one can co-optimize 
the competing objectives of generation cost and security 
indices flEi : see also lIZTl for a stochastic security-constrained 
approach. Albeit computationally complex, stochastic bilevel 
programs are attractive because they can account for the 
coupling between DA and RT (spot) markets ll28l . 1291 . 

All existing MC approaches, however, are centralized. 
Moreover, they are not tailored to address the challenges of 
emerging large-scale dispatchable loads. Specifically, demand 
offers come from demand response (DR) aggregators serving 
large numbers of residential appliances that feature diverse 


Cost function of generator i. 

Utility function of appliance s. 

CVaR transaction cost. 

Sample mean of Fp{-). 

Partial Lagrangian function of the 
stochastic market clearing problem. 
Revenues or payments of the supplier, 
the aggregator, and the wind farm 
located at bus n. 
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Utility functions and inter-temporal constraints. In this context, 
the present paper deals with the DC-OPF based MC with high- 
penetration wind power. Instead of the worst-case or chance- 
constrained formulations, a novel stochastic optimization ap¬ 
proach is proposed to maintain the nodal power balance while 
minimize (maximize) the grid-wide social cost (welfare). The 
social cost accounts for the conventional generation costs, the 
dis-utility of dispatchable loads, as well as a risk measure 
of the cost incurred by (over-) under-estimating the actual 
wind generation. This is essentially a cost of re-dispatching 
the system to compensate wind forecast errors, and is referred 
as transaction cost throughout this paper. The transaction 
cost in the spot market is modulated through an efficient 
risk measure, namely the conditional value-at-risk (CVaR) 
(Sec. IIII-Ab . which accounts not only for the expected cost of 
the recourse actions, but also for their “risks”. A distribution- 
free sample average approximation (SAA) is employed to 
bypass the prohibitively burdensome integration involved in 
the CVaR-based convex minimization ('Sec. lIII-Bli . To clear the 
market in a distributed fashion, a fast and provably convergent 
solver is developed using the ADMM (Sec. llVb . Numerical 
tests are performed to corroborate the effectiveness of the 
novel model and proposed approaches using real power market 
data (Sec. lYli. 

The main contribution of this paper is three-fold: i) a CVaR- 
based transaction cost is introduced for the day-ahead MC 
to judiciously control the risk of (over-) under-estimating the 
wind power generation; ii) a sufficient condition pertinent to 
transaction prices is established to effect convexity of the 
CVaR-based cost; and iii) a distributed solver of the resulting 
stochastic MC task is developed to be run by the market 
operator and DR aggregators while respecting the privacy of 
end users. 

Notation. Boldface lower (upper) case letters represent 
column vectors (matrices); calligraphic letters stand for sets. 

]gd, stand for real spaces of di x d2 matrices, 

d X \ vectors, and non-negative real numbers, respectively; 
Symbols a' and a • b denote the transpose of a, and the 
inner product of a and b; \I\ is the lower endpoint of the 
interval set X. Operator [a]’*' := max{a, 0} is the projection 
to the nonnegative reals, while A (^) indicates the entry-wise 
inequality. Finally, the expectation is denoted by E[-]. 

II. CVaR revisited: A Convex Risk Measure 

Value-at-risk (VaR) and conditional value-at-risk (CVaR) are 
widely used in various real-world applications, especially in 
the finance area, as the popular tools to evaluate the credit risk 
of a portfolio, and reduce the probability of large losses |[30l- 
ll^ . The following revisit is useful to grasp their role in the 
present context. 

Consider a loss function T(x, ^) : AT x S i-A K denoting 
the real-valued cost associated with the decision variable 
X G X C M"; and the random vector ^ with probability 
density function p(^) supported on a set S C K'^. In the 
context of power grids, x can represent the power schedules 
of generators, while ^ collects the sources of uncertainty due 
to for instance renewable energy and forecasted load demand. 



Loss 


Fig. i. Illustration of VaR and CVaR: /3-VaR is the threshold that the loss 
exceeds with at most a small probability 1 — /3. /3-CVaR is the conditional 
expectation of the loss beyond the /9-VaR. 


Clearly, the probability of L(x,^) not exceeding a threshold 
p is given by the right-continuous cumulative distribution 
function (CDF) 

^(x,r;)= J p(Od|. (1) 

Definition 1 (VaR). Given a prescribed confidence level (3 G 
(0,1), the P-VaR is the generalized inverse of 'k defined as 

rip{x) := minlp C M | 'k(x, p) > /?}. (2) 


/3-VaR is essentially the /3-quantile of the random L{x,(f). 
Since T* is non-decreasing in p, (x) comes out as the lower 
endpoint of the solution interval satisfying ik (x, rf) = /3, 
and the commonly chosen values of /3 are, e.g., 0.99, 0.95, 
and 0.9. Clearly, VaR determines a maximum tolerable loss 
of an investment, i.e., a threshold the loss will not exceed 
with a high probability /3. Hence, given the confidence level 
/3, investors are motivated to solve the so-termed portfolio 
optimization problem which yields the optimal investment 
decisions minimizing the VaR value, (x) is proportional to 
the standard deviation if Ti is Gaussian. However, for general 
distributions, /3-VaR is non-subadditive which means the VaR 
of a combined portfolio can be larger than the sum of the 
VaRs of each component. This violates the common principle 
“diversification reduces risk”. Moreover, it is generally non- 
convex rendering the optimization task hard to tackle. 

Because of these conceptual and practical drawbacks, CVaR 
(a.k.a. “tail VaR”, “mean shortfall”, or “mean excess loss”) 
was proposed as an alternative risk metric that has many 
superior properties over VaR. 


Definition 2 (CVaR). The P-CVaR is the mean of the P-tail 
distribution of L{x,^), which is given as 


^'/3(x,?7) := 


0 , 

1-/3 ’ 


ifV< ??/3(x) 
ifV> »7/3(x) 


(3) 
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Truncated and re-scaled from T*, function p is non¬ 
decreasing, right-continuous, and in fact a distribution func¬ 
tion. If Ti is continuous everywhere (without jumps), /3-CVaR 
coincides with the lower CVaR := Eg \L\L > 77 ^(x)], 

that is the conditional expectation of the loss beyond the /3- 
VaR. Hence, roughly speaking, /3-CVaR is the expected loss 
in the worst 100(1 —/?)% scenarios; i.e., cases of such severe 
losses occur only 100(1 — /3) percent of the time. 

The /3-CVaR can be also defined as the optimal value of the 
following optimization problem 

(/)^(x) :=mm|77-f ^-^E^[L(x,|)-? 7 ] + |. (4) 

Let Tjg(x, p) denote the objective function in (|4]i. Key prop¬ 
erties of Fp and its relationship with rjpijx.) and (^/ 3 (x) are 

summarized next. 

Theorem 1 ( OTI . pp. 1454-1457). Function Fp(x, rf) is finite 
and convex in rj. Values P/ 3 (x) and (/>/ 3 (x) are linked through 
Fp{:x.,ri) as 

77/3 (x) = [argmin Fp{x,T])\ (5) 

rjSR 

<(>/3(x) = F/3(x,77/3(x)) ( 6 ) 

min(5fl(x)= min Fb(:x.,ii). (7) 

xex (x,j))gXxR 

Moreover, if L(x,^) is convex in x, then Fp(x,r]) is jointly 
convex in (x, 77 ), while is convex in x. 

From Definition |2] it can be seen that CVaR is an upper 
bound of VaR, implying that portfolios with small CVaR also 
have small VaR. As a consequence of Theorem [T] minimizing 
the convex amounts to minimizing Fp{x,r]), which is 

not only convex, but also easier to approximate. A readily 
implementable approximation of the expectation function Fp 
is its empirical estimate using Ng Monte Carlo samples 
namely 

1 

F/3(x, 77) = 77 -f-^^-^^--^^[L( x,^,)-77]+. ( 8 ) 

Clearly, the sample average approximation method is dis¬ 
tribution free, and the law of large numbers ensures Fp 
approximates well Fp for Ns large enough. Furthermore, 
Fp{'x, Tj) is convex with respect to (x, 77 ) if L(x, is convex 
in X. The non-differentiability due to the projection operator 
can be readily overcome by leveraging the epigraph form of 
F, which will be shown explicitly in Section ITlI-CI 

With the function Fp{x,T]), it is now possible to develop 
the CVaR-based stochastic market clearing, as detailed in the 
next section. 

III. Stochastic Market Clearing 

In a day-ahead electricity market, participants including 
power generation companies and load service entities (LSEs) 
first submit their hourly supply bids and demand offers to 
market operators for the next operating day. Then, the ISO 
or regional transmission organization (RTO) clear the for¬ 
ward markets yielding least-cost unit commitment decisions. 


power dispatch outputs, and the corresponding DA clearing 
prices. The MC procedure proceeds in two stages. A security- 
constrained unit commitment (SCUC) is performed first by 
solving a large-scale mixed integer program to commit gen¬ 
eration resources after simplifying or omitting transmission 
constraints. The second stage involves security-constrained 
economic dispatch (SCED) obtaining the economical power 
generation outputs and the locational marginal prices (LMPs) 
as a byproduct. With unit commitment decisions fixed, SCED 
is usually in the form of DC-OPE, including the transmission 
network constraints 1 ^ . 

The MC process is implemented with a goal of mini¬ 
mizing the system net cost, or equivalently maximizing the 
social welfare. With the trend of increasing penetration of 
renewables, WPPs are able to directly bid in the forward 
market ll^ . Under uncertainty of wind generation, it now 
becomes challenging but imperative for the ISOs/RTOs and 
market participants to extract forecast information and make 
efficient decisions, including reserve requirements, day-ahead 
scheduling, market clearing, reliability commitments, as well 
as the real-time dispatch llTSl . In this section, a stochastic 
MC approach using the CVaR-based transaction cost will be 
developed as follows. 

A. CVaR-based Energy Transaction Cost 

Consider a power system comprising TVf, buses, Ni lines, Ng 
conventional generators, Nyj wind farms and Na aggregators, 
each serving a large number of residential end-users with 
controllable smart appliances. Let T := { 1 , 2 ,..., T} denote 
the scheduling horizon of interest, e.g., one day ahead. If 
a wind farm is located at bus m, two quantities will be 
associated with it: the actual wind power generation Wm, 
and the power scheduled to be injected pw^- Note that the 
former is random, whereas the latter is a decision variable. Eor 
notational simplicity, define also two Wy,-dimensional vectors 
w* := K,..., and := ... ,Pw^J- 

Since w* varies randomly, either energy surplus or shortage 
should be included to satisfy the nodal balance with the 
committed quantity p(y. When surplus occurs, the wind farms 
can sell the excess wind energy back to the spot market, 
or simply curtail it. Eor the case of shortage, in order to 
accomplish the promised bid in the DA contract, farms can 
buy the energy shortfall from the RT market in the form of 
ancillary services. 

Let b* := \b\,... and s* := [s*,..., collect 

the purchase and selling prices at time t, respectively. Clearly, 
with the power shortfall and surplus being [p^.^, — and 
[w* — at time t, the grid-wide net transaction cost is 

T 

r(pw,w) = X] ■ [PVL - w*]+ - s‘ • [w* - p*w]+) 

T 

= • |p‘m/ - W* I -b • (p*^y - w‘)) (9) 

where vj* := and i?* := ; pw and w collect p(y 

and w* for all t € F, respectively. 
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Replacing L{-,-) in (HI with function Fp can be 

expressed through the conditional expected transaction cost as 


F0{pw,r]) =V+Y:7p 


Ev 


i:( 


• \Pw - 


+ • (pV - w*)) -p 


( 10 ) 


A condition guaranteeing convexity of Fp {pw, p) is estab¬ 
lished next. 

Proposition 1. If the selling price does not exceed the 
purchase price for any m € ATy, and t € F, function 
FjsipwjV) jointly convex with respect to (jpw,'n)- 

Proof: Thanks to Theorem [T] it suffices to show that 

T{pw,^) = ELi ■ \Pw - w‘| + - w‘)) 

is convex in pw under the proposition’s condition. Clearly, 
the stated condition is equivalent to ro* ^ 0 for all t € T- 
Thus, by the convexity of the absolute value function, and 
the convexity-preserving operators of summation and expec¬ 
tation ll^ Sec. 3.2], the claim follows readily. ■ 

In this paper, a perfectly competitive market is assumed 
such that all participants act as price takers. That is, every 
competitor is atomistic to have small enough market share so 
that there is no market power affecting the price BtII . For 
American electricity markets, a single pricing mechanism is 
used such that s* = b* holds in most of the scenarios. This 
is a special case of the pricing condition in Prop. [T] which 
facilitates calculating the function (doll since the absolute 
value functions vanish. Note that it is possible that different 
WPPs may buy (sell) wind energy from (to) different sellers 
(purchasers) in a competitive electricity pool as an ancillary 
service, which can yield different purchase and selling prices. 

For most of the European markets including UK, France, 
Italy, and Netherlands, the imbalance prices {b‘,s*}t are 
commonly set in an ex-post way that is known as dual im¬ 
balance pricing ll3^ . Specifically, if the system RT imbalance 
is negative, i.e., the overall market is short, then s‘ = A* ^ b‘ 
holds, where := [Xii • ■ •) collects the DA prices 

at the buses attached with all wind farms. In this case, 
the RT purchase price is typically higher than the DA price, 
reflecting the cost of acquiring the balancing energy 1^ . 
Wind farms with excess energy can sell this part to reduce 
the system imbalance but only be paid the DA prices. On 
the other hand, we have = b* if the market is 

long. Hence, market participants selling excess energy receive 
a balancing price which is lower than the DA one, while 
those running negative imbalance pay the DA price. Note 
that the relationship s* d d b* always holds even when 
the market imbalance outcome is unknown at the time of the 
DA bids. Such a pricing mechanism drives bidders to match 
their forward offers with the true forecasts of generation or 
consumption. 

Leveraging the CVaR-based transaction cost, a stochastic 
MC problem based on the DC-OPF will be formulated next. 


B. CVaR-based Market Clearing 

Let := and — 

[-^DRAi 5 ■ • ■! ^drAn denote the power outputs of the 
thermal generators, and the power consumption of the 
aggregators at slot t, respectively. Define further the sets 
Ma := {l,2,...,Na} and ATg := {1,2,Ng}. Each 
aggregator j G Afa serves a set TZj of residential users, and 
each user r G TZj has a set Srj of controllable appliances. 
Let pjrs be the power consumption of appliance s with 
user r corresponding to aggregator j across the slots. The 
operational constraints of pjrs are captured by a set Vjrs, 
while the end user satisfaction is modeled by a concave 
utility function Ujrs{Pjrs)- Lurthermore, let convex functions 
denote the generation costs, and pgL the base load 
demand. Lor brevity, let vector po collect variables p and 
{PG’PDRA’PW>^*}ter; and vector {Pj}jeMa the power 
consumption of all appliances with the aggregator j. 

Hinging on three assumptions: al) lossless lines, a2) small 
voltage phase differences, and a3) approximated one p.u. 
voltage magnitudes, the DC-OPL based stochastic MC stands 
with the goal of minimizing the social cost: 


T N„ 


Na 


inin^^C'i(PE) - X! UjrsijPjrs') H" 

t—1 2=1 j—lrGTZj, 

S^Sjr 

(11a) 

subject to : 

-^sPg + ^wPw ~ -^oPdra ~ Pel ~ t G F (Hb) 


GAA, f G T 


^g'“ < Pg. < ^gE’ * 

- < PE - Pg~^ ^ ieAfg,tGF 

fmin ^ ^ piax^ ^ ^ j- 

0l = O,t€F 

Odpwd Pw" 

0 < -PdRAj ^ -^DRAj ! j ^ t GF 
^DRA, = jCMa,tCF 

Pjrs C Fjrs-, ty G ^jry ^ C F.j, j G A/a 
variables : {Pj}Eo 


(11c) 

(lid) 

(lie) 

(l lf) 

(llg) 

(llh) 

(Hi) 

(llj) 


where the nodal susceptance matrix B„ := —A^BsA„ G 
jjjg angle-to-flow matrix B/ := —B^A^ G 
row of the branch-node incidence matrix 
An G RMxAffc 2 and —1 in its entry corresponding to 
the from and to nodes of branch i, and 0 elsewhere; and the 
square diagonal matrix Bg := diag(6i,..., is the branch 
susceptance matrix collecting the primitive susceptance across 
all branches. 

Matrices A^ G A^, G and Aa G 

in (fTTbli are the incidence matrices of the con¬ 
ventional generators, the wind farms, and the aggregators, 
respectively. Take Ag as an example, (A^j^n = 1 if the 
nth generator is injected to the mth bus, and {Ag)mn = 0, 
otherwise. Matrices A^, and Aa can be constructed likewise. 
Consider the power network in Pig. |2] adapted from the West¬ 
ern Electricity Coordinating Council (WECC) system ioi. 
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Fig. 2. Modified WECC system featuring 6 buses, 3 generators, 4 aggrega¬ 
tors, 4 base loads, and 3 wind farms. 


With Nf, = 6, Ni = 6, Ng = 3, and Na = 4, matrices Ag, 
Auj, and take the following form: 


'l 

0 

o' 


'l 

0 

o' 


'0 

0 

0 

o' 

0 

1 

0 


0 

1 

0 


0 

0 

0 

0 

0 

0 

1 

A — 

1 -^W — 

0 

0 

0 

5 -^a — 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

0 

0 

0 

0 

0 


0 

0 

1 


0 

0 

1 

0 

0 

0 

0 


0 

0 

0 


0 

0 

0 

1 


A smart appliance example is charging a plug-in hybrid 
electric vehicle (PHEV), which typically amounts to consum¬ 
ing a prescribed total energy Ejrs over a specific horizon from 
a start time to a termination time TJ"/. The consumption 
must remain within a range between p™” and per period. 
With := ..., set Vjrs takes the form: 


^jrs — \ Pjr 


E „t — p. rf a In™™ 

yjrs ^jrs') Fjrs ^ LF'jrs ^Fjrs J5 




V f e p ,, = 0 , Vfer\7't. • (12) 


Further examples of Vjrs and Ujrs {Pjrs 


) can be found in ini, 
is a convex set for several 


where it is argued that Vjrs 
appliance types of interest. 

Linear equality (II Ibb is the nodal balance constraint; i.e., 
the load balance at bus levels dictated by the law of conserva¬ 
tion of power. Limits of generator outputs and ramping rates 
are specified in constraints (II Icb and (IIIdb . Network power 
flow constraints are accounted for in (II leb . Without loss of 
generality, the first bus can be set as the reference bus with zero 
phase in (IIIfb . Constraints (IIIhb and ( 1 lg| i capture the lower 
and upper limits of the energy consumed by the aggregators 
and the committed wind power, respectively. Equality (fTTIli 
amounts to the aggregator-user power balance equation; and 
constraints ( |1 lj| i define the feasible set of appliances. Finally, 
the pre-determined risk-aversion parameter p > 0 controls the 
trade off between the transaction cost and the generation cost 
as well as the end-user utility. 


Remark 1. (Availability of real-time prices). In this paper, 
the real-time prices {b*,s*}tg 7 - are assumed to be perfectly 


known to the ISO for the DA market clearing. However, such 
an assumption can be readily extended to a more practical 
setup by taking the price stochasticity into account. Specifi¬ 
cally, imperfect price information can be modeled by appropri¬ 
ately designing the function T{pw, w) [cf. ^]. For example, 
the expectation can be also taken over the random RT prices 
in (doll as F^(pw,p) = p + T4gE|w.{b‘,s‘}t} “ 

rff'. The dependence between {b*,s*} and w can be further 
investigated. In addition, worst-case analysis is available 
upon postulating an uncertainty set A for {b‘,s‘}. This 
results in a novel risk measure given as L)3(pw,p) = P + 

j4;g]Ew[sup{bt,s*}tGA^(Ptv,w) - p]+. 

It is worth mentioning that SCED and SCUC yield two 
different market pricing systems: locational marginal pricing 
and convex hull pricing (a.k.a. extended LMP). The ED 
formulation produces the LMPs given by the dual variables 
associated with the supply-demand balance constraint. Prices 
supporting the equilibrium solution are found at the intersec¬ 
tion of the supply marginal cost curve with the demand bids. 
However, if discrete operations of UC are involved, there is no 
exact price that supports such an economic equilibrium. This 
issue prompted the introduction of the convex hull pricing 
to reduce the uplift payments HTI . In the present paper, the 
core ED model is considered to deal with the high penetration 
of renewables and large-scale DR programs. Therefore, the 
formulation (fTTb relies on re-solving the dispatch problem with 
fixed UC decisions. 

Remark 2. (Reliability assessment commitment). The pro¬ 
posed dispatch model can be cast as a two-stage program. 
The first stage is the DA MC, and the second is simply the 
balancing operation (recourse action) dealing with differences 
between the pre-dispatch amount and the actual wind power 
generation. Between the DA and RT markets, ISOs implement 
the reliability assessment commitment (RAC) as a reliability 
backstop tool to ensure sufficient resources are available and 
cover the adjusted forecast load online. One principle of the 
RAC process is to commit the capacity deemed necessary to 
reliably operate the grid at the least commitment cost. In this 
step, based on the updated information of the wind power 
forecast, WPPs have an opportunity to feedback to the ISO if 
they are able to commit the scheduled wind power decided by 
the DA MC. Then, the ISO is able to adjust UC decisions as 
necessary to ensure reliability. 

To this end, reformulation of problem (fTTb as a smooth con¬ 
vex minimization is useful for developing distributed solvers, 
as detailed next. 


C. Smooth Convex Minimization Reformulation 

It is clear that under the condition of Proposition [T] the 
objective and the constraints of (fTTb are convex, which renders 
it not hard to solve in principle. Nevertheless, due to the high¬ 
dimensional integration present in Fpi^pwitf) [cf. (fTObl. an 
analytical solution is typically impossible. To this end, it is 
necessary to re-write the resulting problem in a form suitable 
for off-the-shelf solvers. 
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First, as shown in ([8]), an efficient approximation of 
Fpipw^'n) is offered by the empirical expectation using i.i.d. 
samples that is, 


1 r "T 

S—l t —1 


+ ■ [Pw - ^D) - V 


(13) 


Next, by introducing auxiliary variables {us}^Ji, the non¬ 
smooth convex program (HUl can be equivalently re-written as 
the following smooth convex minimization; 

T Ng Ng, 

min^^C'i(P^J Y UjrsiPjTs) 
t —1 i —1 j—lr^TZj, 

S^Sjr 

+ '■(’' +fey) 

subject to : (fTTbli — 

T 

Y (^* ■ I -f # ■ {Pw - w‘)) <Us+v, 

t^l 

s€j\fs (14b) 

variables : {Pj}Y, {us S 

Under mild conditions, the optimal solution set of (fT4l i con¬ 
verges exponentially fast to its counterpart of (fTTTi . as the 
sample size Ng increases. The proof is based on the theory 
of large deviations ll42l . but is omitted here due to space 
limitations. 

Problem (fT4l) can be solved centrally at the ISO in principle. 
However, with large-scale DR, distributed solvers are well 
motivated not only for computational efficiency but also for 
privacy reasons. Specifically, functions Ujrs{Pjrs) and sets 
{'Pjrs} are private, and are not revealed to the ISO; and (ii) 
the operational sets {'Pjrs}j,r,s of very large numbers of het¬ 
erogenous appliances may become prohibitively complicated; 
e.g., mix-integer constraints can even be involved to model the 
ON/OFF status and un-interruptible operating time of end- 
user appliances ll43l . Il44ll . This renders the overall problem 
intractable for the ISO. To this end, the DR aggregators can 
play a critical role to split the resulting optimization task as 
detailed next. 


IV. Distributed Market Clearing via ADMM 

Selecting how to decompose the optimization task as well 
as updating the associated multipliers are crucial for the 


distributed design. Fewer updates simply imply lower com¬ 
munication overhead between the ISO and the aggregators. 
One splitting approach is the dual decomposition with which 
the dual subgradient ascent algorithm is typically very slow. 
Instead, a fast-convergent solver via the ADMM ll45ll is 
adapted in this section for the distributed MC. 


A. The ADMM Method 

Consider the following separable convex minimization prob¬ 
lem with linear equality constraints: 


min /(x)-f g(y) ( 15 a) 

subject to : Ax -f By = c. ( 15 b) 


For the stochastic MC problem (0, the primal variable 
X comprises the group {msIsgVs po, while y collects 
{Pj}j^Mg- Hence, set X captures constraints (IIlbl i- (llIhl) 
and (I14bl i while y represent s (|llj| l. The linear equality con¬ 
straint (I15bb corresponds to (IIlib . 

Let A := [A},..., denote the Lagrange 

multiplier vector associated with the constraint (II lib . The 
partially augmented Lagrangian of (fl4b is thus given by (fThb . 
where the weight p > 0 is a penalty parameter controlling 
the violation of primal feasibility, which turns out to be the 
step size of the dual update. As the iterative solver of (fThb 
proceeds, the primal residual converges to zero that ensures 
optimality. Judiciously selecting p thus strikes a desirable 
tradeoff between the size of primal vis-a-vis dual residuals. 
Note also that by varying p over a finite number of iterations 
may improve convergence ll45l . In a nutshell, finding the 
“optimal” value of p is generally application-dependent that 
requires a trial-and-error tuning. 

Different from Bbl where the power balance and phase 
consistency constraints are relaxed, in this work only the 
aggregator-user power balance equation (fTTIb is dualized so 
that the nodal balance equation (II Ibb is kept in the subproblem 
of the ISO. Decomposing the problem (fT4b in such a way 
can reduce the heavy computational burden at the ISO while 
respect the privacy of end users within each aggregator. 
The ADMM iteration cycles between primal variable updates 
using block coordinate descent (a.k.a. Gauss-Seidel), and dual 
variable updates via gradient ascent. The resulting distributed 
MC is tabulated as Algorithm[Tl where k is the iteration index. 
The last step is a reasonable termination criterion based on the 
primal residual ll45l Sec. 3.3.1] 


e:= 


T Na / 
t=lj=l \ 



(17) 


T Ng Ng / ^Ns 

y, A) = ^ ^ ^ Y^ Ujrs{.P]rs) AT Q) 

j=lreTlj, V *'■ 

S^Sjr 


t^l 2=1 


T Na ( 

' ^ ^ ( -^dra,- - Y P 

t —1 j=l \ r,s 


* 1 -I- - 

'jrs I “E 2 


T Ng 


9 X) X) ( ^DRAj -Yp 

t=l i=l \ 


'jrs 


r.s 


r,s 


(16) 
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Algorithm 1 ADMM-based Distributed Market Clearing 

1: 

Initialize A(0) = 

0 


2: 

repeat for k = 

0,1,2,. 


3: 

update primal variables: 


x(fc + 1) = 

arg min 

-(^p(x,y(fc), A(fc)) (18) 



xGA' 



y(fc + 1) = 

arg min 

Lp(x(fc -b l),y, A(fc)) (19) 



yey 


4: 

update dual 

variables: for all j G Ma and t eT 


(fc -P 1) _ A‘(/c) -I- p(PDRAj (^ + 1) ~ Y.P^rsjf^ + 1)) 




(20) 

5: 

until ^ 




Specifically, given the Lagrangian multipliers A(fe) and the 
power consumption {pjrs(^)}jr-s of the end-user appliances, 
The ISO solves the convex subproblem (fTSl l given as follows: 


convex problem (l2^ . each aggregator must collect the corre¬ 
sponding users’ information including Ujrs and Vjrs- This is 
implementable via the advanced metering infrastructure El. 

Remark 3. (Distributed demand response). It must be further 
pointed out that the quadratic penalty (TdrAj ~ 'Yhr sP^jrs) 
in (O couples load consumptions {p*>s} over different resi¬ 
dential users. Hence, the ADMM-based distributed solver may 
not be applicable whenever must be updated per end 
user rather than the aggregator. This may arise either to 
strictly protect the privacy of end users from DR aggregators, 
or, to accommodate large-scale DR programs where each 
aggregator cannot even afford solving the subproblem (Hal). 
In this case, leveraging the plain Lagrangian function (no 
coupling term), the dual decomposition based schemes can be 
utilized by end users to separately update {p^jrs} parallel; 
see e.g., snu and l\48\l . 

The convergence of the ADMM solver and its implications 
for the market price are discussed next. 


Po(fc + l) = argmin ^ +p(r]+ ) 

po.{«a ttr. V 

+ ^ •^5(^)^DRAj ^ f-PoRAj ~ I 

teT, teT, \ r,s ) 

(21a) 

subject to : 

-^sPg ^wPw ~ ■'^oPdra ~ Pel “ t G P (21b) 

Pg-"" < Pg, < Paf’ * e A4, ter ( 21 c) 

- Rf™ < - P^-^ < R"P, ieMg,teT (2ld) 

fmin ^ ^ finax^ ^ ^ ^ ( 21 e) 

0 * = 0 , f G r ( 2 if) 

Oy(pwe< Pw" (21g) 

0 < PuKAj ^ -PdRAjI j S t eT (21h) 

T 

■ \p\^ - w‘| -b • (p*^^ - w‘)) <Us+V, 

t=l 

and Us > 0, s G A/’s- (21i) 

Interestingly, ( fT9b is decomposable so that {Pj>s(/^)}r,s can 
be separately solved by each aggregator: 


{Pj>s(A: + 1 )}as = argmin - Y ^^jWYPjrs 

{Pjrs}r,s 


t —1 

T / 


“ Y^ ^psiPjrs) + ^ X/ ( Y^Pps - PoRAj (k + 1) 


s^Sjt 


t —1 \ r.s 


subject to . {pjrs € ^jrs}r,s- 


(22a) 

(22b) 


Having found Po(fc) and {Pj>s(/i^)}jrs> the multipliers {Pj}j^t 
are updated using gradient ascent as in (l20l) . To solve the 


B. Pricing Impacts 

Suppose two additional conditions hold for the convex 
problem (fT4t : cl) functions {Ci{-)}i and {—Ujrs{-)}jrs are 
closed and proper convex; and c2) the plain Lagrangian Lq has 
a saddle point. Then, the ADMM iterates of the objective (I14al) 
and the dual variables are guaranteed to converge to the 

optimum ll45l . In addition, if the objective is strongly convex, 
then the primal variable iterates including pc, Pdra, Pw and 
{Pj}jGHa. converge to the globally optimal solutions. 

The guaranteed convergence of the dual variables also facil¬ 
itates the calculation of LMPs. Let A* := [A*,..., ]' and 

T* := [fl,... denote the optimal Lagrange multipliers 

associated with the aggregator-user balance constraint (II lit , 
and the nodal balance constraint (II Ibk respectively. Note that 
with the optimal solutions A* and {pjrsjjVs obtained by the 
ADMM solver, the LMPs {t*}^ can be found by solving the 
subproblem (l2Tl i with primal-dual algorithms. In addition, if 
0 PdrAj ^ -^DRAj! t holds at the optimal solution 
Pqra , then A* = AJjT*; i.e.. A* = f* for all aggregators 
j attached with bus n (see also iflTl ). To this end, payments 
of the market participants can be calculated with the obtained 
LMPs and optimal DA dispatches. In the RT market of a two- 
settlement system, if the supplier at bus n delivers Pq^ with 
the real-time price f^, then the supplier gets paid 

T 

nG„ = E 

Likewise, the aggregator at bus n needs to pay 

T 

HdRA^ = ^ ^n^DRAn + ^ni^DRAr^ ~ ^DRA„)- 
t^l 

The revenue of the wind farm at bus n is 

T 

niV„ = Y i^nPw„ + Snbvl-pwY “ K[Pw„ “ 
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TABLE I 

Conventional generator parameters. The units of a; and bi are 
$/(MWH)^ and $/MWH, respectively. THE REST ARE IN MW. 


Unit 

(Li 

h 

pmax 

^Gi 

pmin 

^Gi 

Rr 

pdown 

1 

0.3 

50 

90 

10 

50 

50 

2 

0.15 

30 

50 

5 

35 

40 

3 

0.2 

40 

60 

8 

40 

40 


TABLE II 

Parameters of PHEVs. All listed hours are the ending ones; 

W.P. MEANS WITH PROBABILITY. 


SpHBV (kWh) 

Uniform on {10, 11, 12} 

PpSev (kWh) 

Uniform on {2.1, 2.3, 2.5} 

PPHBV (kWh) 

0 

T^St 

jrs 

lam 

T-iend 

jrs 

6am w.p. 70%, 7am w.p. 30% 


Remark 4. (Pricing consistence). In a perfectly competitive 
market, any arbitrage opportunities between the DA and RT 
markets are exploited by market participants. Hence, the DA 
nodal prices are consistent with the DT nodal prices meaning 
the expectations of the latter converge to the former. The 
concepts of price distortions and revenue adequacy have been 
recently proposed for the stochastic MC in H49\l . In the setup 
of a single snapshot therein, it has been proved that the 
medians and expectations of RT prices converge to the DA 
counterparts for the and li penalties between the RT and 
DA power schedules, respectively. Building upon this solid 
result, it is possible to establish bounded price distortions for 
the proposed model, while its consistent pricing property can 
also be analyzed in a similar fashion. The involved important 
analysis is however beyond the scope of this paper, and is left 
for future work. 

V. Numerical Tests 

In this section, simulated tests are presented to verify the 
merits of the proposed CVaR-based MC. The tested power 
system is modified from the WECC system as illustrated in 
Fig. |2] Each of the 4 DR aggregators serves 200 residential 
customers. The scheduling horizon starts from 12am until 
23pm, a total of 24 hours 

Time-invariant generation cost functions were chosen 
quadratic as Ci{PQ.) = ai{PQ.)‘^ + biPQ. for all i and 
t. For simplicity, each end user has one PHEV to charge 
from midnight. All detailed parameters of the conven¬ 
tional generators and loads are listed in Tables U and |II] 
The upper bound of each aggregator’s consumption is 
PbrAj = 50 MW. At a base of 100 MVA, the values of 
the network reactances are {Xie, X 62 , -A25, ^53, ^34, X41} = 
{0.2,0.3,0.25,0.1,0.3,0.4} p.u. Finally, no flow limits were 
imposed, while the utility functions {C{jrs(')} were set to 
zero. The resulting convex programs (EB and (I 22 I 1 were 
modeled using the Matlab-based package CVX ESI, and solved 
by SeDuMi lISTl . 

Variable characteristics of the daily power market are cap¬ 
tured via two groups of parameters shown in Fig. |3] the fixed 
base load demand {PblI’ ^nd the purchase prices {b*} at 
the buses attached with three wind farms. The prices were 




Fig. 3. Fixed base load demand {p|l} energy purchase prices {b*}. 



Fig. 4. Convergence of the objective value (Hi and the primal residual {nj. 

obtained by scaling the real data from the Midcontinent ISO 
(MISO) II 52 I . Two peaks of {b*} appear during the morning 
7am to 12pm, and early night 6pm to 9pm. The selling prices 
{s*} were set to s* = 0.9b‘ satisfying the convexity condition 
in Proposition [T] The rated capacity of each wind farm was 
set to 20 MW, yielding a 23% wind power penetration of the 
total power generation capacity. 

Wind power output samples {w}}s t are needed as inputs 
oflEB. These samples can be obtained either from forecasts of 
wind power generation, or, by using the distributions of wind 
speed together with the wind-speed-to-wind-power mappings 
[cf. El]. In this paper, the needed samples were obtained from 
the model = w* -f n‘, Vf G T. The DA wind power 
forecasts {w*} were taken from the MISO market on March 
8, 2014. The forecast error n} was assumed zero-mean white 
Gaussian. Possible negative-valued elements of the generated 
samples were truncated to zero. Finally, the sample 

size Ng = 200, the probability level /3 = 0.95, the trade-off 
weight p, = 1, and the primal-residual tolerance eP’'' = lO”"* 
were set for all simulations, unless otherwise stated. 

Figure |4] demonstrates the fast convergence of the proposed 
ADMM-based solver. The pertinent parameters were set to 
p = 35 and A* (0) = p*>s(0) = 0. Clearly, both the cost and 
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Optimal total cost {$) 


Fig. 5. Empirical CDFs of the optimal social cost. 

TABLE III 

Mean and standard deviation of the total cost and the 

CONVENTIONAL GENERATION COST: RISK-LIMITING VERSUS NO 
RISK-LIMITING DISPATCH. THE UNITS ARE ALL IN $. 


Dispatch scheme 

Mean 

Std 

Conv. gen. cost 

CVaR-based risk-limiting 

44363.26 

493.15 

26047.66 

With expected wind power 

50095.68 

498.13 

50194.59 

Without wind power 

51619.24 

476.25 

57122.82 


the primal residual converge very fast to the optimum within 
10 iterations. Note that due to the infeasibility of the iterates 
at the beginning, the objective function starts from a value 
smaller than the optimum, and then monotonically converge 
to the latter. 

Three methods were tested to show the performance of 
the optimal dispatch and cost; (i) the novel CVaR-based risk- 
limiting MC; (ii) the no risk-limiting MC with the expected 
wind power generation {w*}; and (iii) the MC without wind 
power integration. Specifically, = w* was simply used 
in the nodal balance (I21bl) for (ii), while = 0 for 
(iii). There are no CVaR-pertinent terms in the objective 
and constraints for the last two alternatives. For all three 
approaches, the generation cost J2t=i Si=i ) fixed 

after solving (fT4l i. Hence, randomness of the optimal total cost 
stems from the transaction cost due to the stochasticity of the 
actual wind power generation {w*} [cf. ®]. 

Figure |3 presents the cumulative distribution functions 
(CDFs) of the optimal total costs using 100, 000 i.i.d. wind 
samples with mean {w*}. Clearly, the two competing alter¬ 
natives always incur higher costs than the novel CVaR-based 
approach. The values of the mean and standard deviation (std) 
of the optimal total cost are listed in Table HII] It can be 
seen that, compared with the other two methods, the proposed 
scheme has a markedly reduced expected total cost and small 
changes in the std. 

Figures |6] |7] and |8] compare the optimal power dispatches 
{PG’Pvu> PoRAiiGT of the proposed scheme with those of 
the scheme (ii). In Fig. |6] it can be clearly seen that over a 
single day the CVaR-based MC dispatches lower and smoother 
Pg than the one with (ii). Furthermore, for the novel method, 
generators 1 and 3 are dispatched to output their minimum 
generation P™", while the output of the generator 2 changes 
within its generation limits across time. Such a dispatch results 




12am 3am 6am 9am 12pm 3pm 6pm 9pm 12am 


No risk-limiting dispatch with expected wind power 
Fig. 6. Optimal power dispatch of 




Fig. 7. Optimal power dispatch of pw". 


from the economic incentive since the unit 2 has the lowest 
generation cost among all three generators [cf. Table [U. On the 
contrary, both generators 2 and 3 fluctuate within a relatively 
large range in (ii), mainly to meet the variation of base load 
demand pbl\ see Fig. [3] 

As shown in Figure |7] the novel CVaR-based approach also 
dispatches more p|y than that of (ii). This is because the 
energy purchase prices b* are smaller than the conventional 
generation costs [cf. Table U and Fig. |3. In addition, and 
contribute most of the committed wind power at 1pm and 
2pm due to the cheaper buying prices during the corresponding 
slots [cf. Fig.|3l. Interestingly, Figure [8] shows that the PHEVs 
are scheduled to start charging earlier for the CVaR-based MC, 
where Pdra is jointly optimized with pG and pvu- 

Finally, Figure |9] shows the effect of the weight parameter 
fi on the optimal costs of the conventional generation and 
the CVaR-based transaction. As expected, the CVaR-based 
transaction cost decreases with the increase of p. For a larger 
p, less p|y is scheduled so that more wind power is likely to 
be sold in the RT market that yields selling revenues rather 
than purchase costs. Consequently, to keep the supply-demand 
balance, higher conventional generation cost is incurred by the 
increase of p^. 
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^^■Aggregator-1 
I I Aggregator-2 

I I Aggregator-3 

Aggregator-3 



2am 3am 4am Sam 6am 
CVaR-based risk-limiting dispatch 



Fig. 8. Optimal power dispatch of pdra- 



0.001 1 2 10 




Fig. 9. Optimal costs of conventional generation and CVaR-based transaction. 


VI. Conclusions and Future Work 

Day-ahead stochastic market clearing with high-penetration 
wind power was investigated in this paper. A stochastic 
optimization problem was formulated to minimize the market 
social cost consisting of the generation cost, the utility of 
dispatchable loads, as well as the CVaR-based transaction 
cost. The SAA method was introduced to bypass the inher¬ 
ent high-dimensional integral, while an ADMM-based solver 
was developed to clear the market in a distributed fashion. 
Extensive tests on a modihed WECC system corroborated the 
effectiveness of the novel approach, which offers risk-limiting 
dispatch with considerably reduced conventional generation. 

A number of appealing directions open up towards ex¬ 
tending the proposed framework. Eirst, it is interesting to 
study the extended LMPs by solving a large-scale stochastic 
SCUC with start-up (-down) and no-load costs. Second, a 
deep explore of the price consistence for multi-period time¬ 
coupling MC is in our research agenda. Additional topics 
worth further investigation include congestion management, 
reserve procurement, as well as security assessment issues. 
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